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Abstract 

A  numerical  algorithm  and  the  corresponding  paralleled  implementation  for  the  study  of  magnetohydrodynamics  (MHD)  of 
large  density  ratio,  three-dimensional  multiphase  flows  at  low  magnetic  Reynolds  numbers  have  been  developed.  The  algorithm 
employs  the  method  of  front  tracking  for  the  propagation  of  material  interfaces  and  the  embedded  interface  method  for  solving 
elliptic-parabolic  problems  associated  with  approximations  of  incompressible  fluids  and  low  magnetic  Reynolds  numbers.  The 
use  of  embedded  interface  method  supports  arbitrary  discontinuities  of  density  and  other  physics  properties  across  interfaces 
and  significantly  improves  methods  that  smear  interface  discontinuities  across  several  grid  cells.  The  numerical  algorithm  has 
been  implemented  as  an  MHD  extension  of  FronTier,  a  parallel  front  tracking  hydrodynamic  code,  verified  using  asymptotic 
solutions  and  validated  through  the  comparison  with  experiments  on  liquid  metal  jets.  The  FronTier-MHD  code  has  been  used 
for  simulations  of  liquid  mercury  targets  for  the  proposed  muon  collider  /  neutrino  factory,  ablation  of  pellets  in  tokamaks,  and 
processes  in  hybrid  magnetoinertial  fusion. 

Keywords:  Front  tracking,  multiphase  MHD,  liquid  metal  MHD 


1.  Introduction 

MHD  of  liquid  metals  and  conducting  liquid  salts  attracts  considerable  attention  from  researchers  because  of 
their  current  and  potential  applications  in  fusion  energy  research,  accelerator  sciences,  and  industrial  processes. 
Applications  such  as  liquid  wall  plasma  facing  components  (PFCs)  in  devices  for  magnetically  confined  fusion 
(tokamaks)  [1]  -  [2],  liquid  metal  targets  for  future  particle  accelerators  [3],  etc.,  motivate  the  study  of  free  surface 
magnetohydrodynamic  flows  either  in  vacuum  or  non-conducting  media. 

Among  other  codes  developed  to  study  such  physical  problems,  we  would  like  to  single  out  HIMAG  and  com¬ 
pressible  FronTier-MHD.  The  HIMAG  code  [4]  (HyPerComp  Incompressible  MHD  solver  for  Arbitrary  Geome¬ 
tries)  is  developed  to  model  the  flow  of  liquid  metal  with  free  surfaces  in  the  presence  of  strong  multi-component 
magnetic  fields.  It  uses  hybrid  meshes  comprising  of  hex,  prism  and  tet  elements,  and  the  level  set  technique 
for  the  free  surface  support.  But  the  diffusive  nature  of  the  level  set  method  that  replaces  the  density  discontinu¬ 
ity  across  material  interfaces  by  a  continuous  density  function  and  smears  the  discontinuity  across  several  mesh 
blocks  limits  the  accuracy  for  problems  involving  large  density  discontinuities. 
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The  compressible  FronTier-MHD  code  [5],  developed  by  one  of  the  authors  and  collaborators  as  an  extension 
of  the  hydrodynamic  code  FronTier  [11],  employs  the  method  of  front  tracking  [13]  for  material  interfaces,  second 
order  hyperbolic  solvers  for  equations  of  compressible  fluids,  and  the  embedded  boundary  method  for  the  elliptic 
problem  in  complex  domains.  The  FronTier  code  always  keeps  discontinuities  sharp  and  eliminates  or  strongly 
restricts  numerical  diffusion  across  material  interfaces.  It  supports  large  number  of  geometrically  complex  inter¬ 
faces  in  two-  and  three-dimensional  spaces  and  robustly  resolves  their  topological  changes.  FronTier  has  been 
widely  used  for  variety  of  fundamental  science  (turbulent  fluid  mixing  [6])  and  applied  problems  (liquid  targets 
for  particle  accelerators  [7],  pellet  ablation  in  tokamaks  [8],  and  plasma  jet  liners  for  magneto-inertial  fusion  [9]). 
The  code  is  well  suited  for  free  surface  MHD  phenomena  driven  by  hydro  waves,  for  instance  in  the  case  of  matter 
interacting  with  strong  energy  sources.  The  obvious  limitation  of  this  code  for  the  simulation  of  slow  flows  of 
liquid  metals  is  the  restriction  of  time  steps  by  the  CFL  condition  due  to  acoustic  waves. 

In  order  to  overcome  this  limitation  but  keep  all  advantages  of  front  tracking,  we  have  developed  a  sharp 
interface  MHD  algorithm  for  incompressible  multiphase  MHD  flows  in  the  low  magnetic  Reynolds  number  ap¬ 
proximation.  As  the  method  is  dependent  on  the  quality  of  the  Navier-Stokes  equation  solver,  we  would  like  to 
comment  first  on  the  hydrodynamic  component  of  the  code.  Front  tracking  has  already  been  used  for  the  simula¬ 
tion  of  incompressible  Navier-Stokes  equations  [14,  15].  But  unlike  the  front  tracking  method  for  compressible 
flows  [11]  which  always  keeps  the  density  discontinuity  sharp,  previous  implementations  of  the  front  tracking  for 
incompressible  flows  employed  the  smoothing  of  density  similar  to  the  level  set  method.  Other  methods  such 
as  the  ghost  fluid  method  [16,  17]  and  the  immersed  interface  method  [18,  19]  also  have  difficulties  with  large 
density  ratios  across  the  interface.  A  front  tracking  algorithm  for  incompressible  Navier-Stokes  approximations 
that  successfully  deals  with  the  large  density  discontinuity  problem  has  been  recently  proposed  by  authors  and 
collaborators  [20]  by  using  the  embedded  boundary  method  [21]  . 

In  this  paper,  we  describe  a  front  tracking  MHD  algorithm  for  free  surface  /  multiphase  flows  in  the  low 
magnetic  Reynolds  number  approximation,  coupled  with  the  incompressible  hydrodynamic  solver,  as  well  as 
validation  and  verification  tests.  With  the  advantage  of  incompressible  hydrodynamic  solver,  the  code  can  deal 
with  the  simulation  of  large  time  scales,  in  particular,  with  flows  of  free  surface  liquid  metals  in  magnetic  fields. 

The  paper  is  organized  as  follows.  The  governing  equations  and  approximations  are  discussed  in  Section  2. 
In  Section  3,  the  numerical  algorithm  and  implementation  are  described.  The  verification  and  validation  test  is 
presented  in  Section  4.  We  briefly  describe  applications  of  the  code  in  the  area  of  accelerator  targets  in  Section  5. 
Finally,  we  conclude  the  paper  withe  the  summary  of  our  results  and  perspectives  for  the  future  work. 

2.  Governing  Equations 

We  are  interested  in  the  description  of  multiphase  or  multi-material  systems  involving  conducting  fluids  inter¬ 
acting  with  neutral  fluids  or  gases  in  the  presence  of  magnetic  fields.  Interfaces  of  the  phase  or  material  separation 
are  assumed  to  be  sharp  (the  thickness  of  the  interface  is  negligible)  and,  in  general,  geometrically  complex.  The 
numerical  simulation  of  liquid  metals,  liquid  salts,  and  weakly  ionized  plasmas,  which  are  relatively  weak  electri¬ 
cal  conductors,  is  difficult  using  the  standard  full  systems  of  MHD  equations  [12].  Fast  diffusion  of  the  magnetic 
field,  caused  by  low  value  of  electrical  conductivity,  introduces  unwanted  small  time  scales  into  the  problem.  If 
the  time  scale  of  the  diffusion  of  the  magnetic  field  is  small  compared  to  hydrodynamic  time  scale,  the  magnetic 
Reynolds  number  [24] 

Re«  =  hf?!, 

cz 

where  L  is  the  typical  length  scale,  u  is  the  fluid  velocity,  and  cr  is  the  electric  conductivity,  is  small.  If,  in  addition, 
the  eddy-current-induced  magnetic  field  5B  is  small  compared  to  the  external  field  B,  the  full  system  of  MHD 
equations  can  be  simplified  by  neglecting  the  time  evolution  of  the  magnetic  field.  In  this  case,  the  generalized 
Ohm’s  law  is  used  for  the  evaluation  of  the  current-density  distribution  instead  of  the  Maxwell  equation  J  = 
^VxH,  where  the  magnetic  field  H  and  the  magnetic  induction  B  are  related  by  the  magnetic  permeability 
coefficient  //:  B  =  /./II.  The  governing  equations  of  incompressible  conductive  fluids  in  the  low  magnetic  Reynolds 
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number  approximation  are 


V  u 

J 

V- J 


yuAu  -VP  +pg  +  -(J  X  B) 
c 

0 

cr|-V</3  4 - U  X  fij 

0 


(1) 

(2) 

(3) 

(4) 


Taking  the  divergence  of  both  sides  of  equation  (3)  together  with  equation  (4),  an  elliptic  equation  for  the  electric 
potential  is  obtained: 

V  •  (crVip)  =  V--(uxB)  (5) 

c 

If  a  conductive  fluid  interfaces  a  neutral  fluid  or  gas,  the  current  density  vector  is  tangential  to  the  material  inter¬ 
face.  This  statement  is  expressed  by  the  following  Neumann  boundary  condition  for  the  Poisson  equation  (5). 


dip 
<9n  r 


-(u  x  B)  ■  n 

c  r 


(6) 


where  T  is  the  boundary  of  conductive  fluid. 


3.  Numerical  Algorithm  and  Implementation 

The  proposed  numerical  algorithm  uses  the  method  of  front  tracking  for  the  propagation  of  material  interfaces 
and  the  embedded  interface  method  for  elliptic  problems  associated  with  an  implicit  discretization  of  the  incom¬ 
pressible  Navier-Stokes  equations  (1)  -  (2)  and  the  Poisson  problem  for  the  electric  potential  (5)  -  (6).  The  task  is 
complicated  by  the  fact  that  these  elliptic  problems  contain  either  geometrically  complex  outside  boundary  or  an 
interior  surface  across  which  large  discontinuities  of  material  properties  or  solutions  occur.  Within  the  method  of 
front  tracking,  the  fluid  interface  is  represented  as  an  explicit  co-dimension  one  Lagrangian  mesh  moving  through 
a  volume  filling  Eulerian  mesh. 


3.1.  Elliptic  Interface  Problem 

The  embedded  boundary  method  (EBM)  for  irregular  domains  [21,  22]  was  extended  by  the  authors  ([26]) 
to  solve  the  elliptic  and  parabolic  interface  problem  with  interior  boundaries.  In  general  form,  the  elliptic  partial 
differential  equation  is 

V-/3VT  =  /  (7) 

where  the  function  f>  is  continuous  and  smooth  except  at  the  interior  boundary  and  f  is  a  given  function,  continuous 
except  perhaps  at  the  interior  boundary.  Boundary  conditions  for  both  exterior  and  interior  boundary  are  needed 
to  close  the  problem.  The  boundary  condition  for  the  exterior  boundary  is  either  of  Dirichlet  or  Neumann  type. 
At  the  interior  boundary,  the  following  jump  conditions  are  applied: 


IT] 

dp 


J\ 


(8) 

(9) 


The  embedded  boundary  method  is  a  finite  volume  method  for  an  irregular  domain  embedded  on  a  Cartesian 
grid.  When  embedded  boundary  method  is  used  to  solve  the  elliptic  boundary  value  problem,  unknowns  are 
defined  in  the  computational  cell  centers  for  both  interior  cells  (full  cells)  as  well  as  boundary  cells  (partial  cells) 
which  intersect  with  the  interior  boundary,  instead  of  being  defined  in  the  geometrical  center  of  a  cell  ([26,  21, 
22,  25]).  This  main  feature  of  EBM  is  retained  in  the  algorithm  for  the  elliptic  interface  problem.  For  a  full  cell, 
one  unknown  is  set  at  the  cell  center  and  the  standard  finite  volume  method  is  used  to  obtain  one  linear  algebraic 
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(a) 


(b) 


(c) 


Fig.  1.  (a)  Placement  of  unknowns  in  a  cell  containing  the  interface,  (b)  stencil  for  the  interface  unknowns  for  the  jump  condition,  and  (c) 
stencil  for  the  cell  center  unknowns. 


equation.  On  the  other  hand,  more  unknowns  are  needed  in  order  to  discretize  the  elliptic  equation  consistent  with 
two  interface  jump  conditions  for  partial  cells.  Figure  la  shows  the  placement  of  unknowns  in  a  cell  with  interior 
boundary.  The  whole  cell  contains  two  partial  cells  representing  two  material  components  (A  and  B)  which  are 
separated  by  the  interior  boundary.  Two  unknowns  (7j,  T2)  are  defined  in  the  cell  center  for  each  part  of  the 
cell.  In  order  to  satisfy  the  jump  conditions  (8)  and  (9),  two  additional  unknowns  Ta  and  7*  are  defined  in  the 
geometrical  center  of  the  segment  of  the  interior  boundary  intersecting  with  the  cell. 

While  the  algorithm  implemented  in  the  code  works  in  both  two-  and  three-dimensional  spaces,  we  focus  here 
on  a  two-dimensional  case  for  visual  simplicity.  A  schematic  of  the  corresponding  stencil  for  the  interpolation  of 
boundary  unknowns  is  shown  in  Figure  lb.  The  direction  of  the  normal  vector  to  the  interior  boundary  is  from  A 
to  B,  assuming  that  A  is  the  interior  component.  The  discretization  of  the  jump  condition  (8)  is  simply 

Ta  ~Tb  =  (10) 


The  discretization  of  the  jump  condition  (9)  is  more  complicated  because  the  normal  derivatives  of  unknowns 
in  both  sides  of  the  interior  boundary  have  to  be  calculated.  This  is  processed  by  fitting  a  quadratic  polynomial 
([21,  22]).  The  main  idea  is  to  use  one  unknown  in  cell  (i,j),  two  unknowns  in  the  first  layer  of  neighbors  of  cell 
(i,j),  and  three  unknowns  in  the  second  layer  of  neighbors  of  cell  (i,j),  that  six  unknowns  in  total  are  given  for  six 
coefficients  of  the  quadratic  polynomial.  For  example,  in  order  to  construct  the  quadratic  polynomial  to  evaluate 
the  flux  at  the  interior  boundary  segment  center  for  component  A,  the  six  unknowns  are  Ta,  T,+\  j,  Ti+i,j-\,  Ti+2j, 
Ti+2,j-u  and  Similarly,  we  can  construct  the  quadratic  polynomial  for  component  B.  Taking  the  normal 

derivatives  of  the  fitted  polynomials  for  two  components  to  obtain  §^|A,  respectively,  and  using  the  jump 
condition  (9),  we  obtain. 


P\b 


3T 

dn 


dT 

B  dn 


h 


The  embedded  boundary  method  is  used  to  setup  two  equations  for  two  unknowns  in  the  cell  center 
For  the  partial  cell  cdef,  a  similar  stencil  is  used  to  discretize  the  elliptic  operator  (see  Figure  lc). 
equation  (7)  and  using  the  divergence  theorem,  we  obtain 


(11) 

separately. 

Integrating 


^cdef 


V  ■  (fiVT)ds 


-i 

Jd( 


dicdef) 


pVTndl 


which  is 


n  dl  + 


The  discretized  form  is 


led  •  FlllXcd  +  Ide  ■  FluXde  +  lef  ■  Fluxef  +  lfC  ■  FluxfC 


(12) 


(13) 


ds 


(14) 
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where  lxy  is  the  length  of  the  segment  between  x  and  y.  For  Fluxrd,  a  second  order  derivative  is  calculated  by 


and 


using  the  linear  interpolation  of 
For  Fluxde ,  we  simply  use  central  difference 


Tj+ij-i-Tj+ij 


Ti+i  j-Tt 
Ay 


Ax 


in  the  center  of  segment  cd  [21]  and  multiplying  by  [i. 
-  to  calculate  the  derivative  and  multiply  fl.  Fluxef  is  obtained 


j1t  _  _ ^ 

similarly  to  Fluxcd  by  evaluating  linear  interpolation  of  IJ*^  IJ  and  l+1,-'H^  '+lj  in  the  center  of  ef  and  multiplying 
by  /3.  Fluxfc  is  evaluated,  as  described  in  the  previous  paragraph,  as  fi\A  |^|  .  Similarly,  fluxes  of  other  partial 
cells  are  obtained.  As  the  unknowns  at  the  geometrical  centers  of  interface  segments  can  be  expressed  in  terms  of 
unknowns  in  cell  centers,  the  resulting  system  of  equations  is  written  only  for  unknowns  defined  at  cell  centers. 

The  3D  algorithm  is  similar.  A  bi-linear  interpolation  is  used  for  the  interface  fluxes  and  10  unknowns  are 
used  to  construct  a  quadratic  polynomial  in  3D  containing  10  coefficients. 


3.2.  Flydro-  and  MFID  Algorithms 

The  system  of  MHD  equations  ( 1 )  -  (4),  a  coupled  parabolicelliptic  system  in  a  geometrically  complex  domain, 
is  solved  using  operator  splitting  and  front  tracking.  The  propagation  and  redistribution  of  the  interface  using  the 
method  of  front  tracking  ([11], [13])  is  performed  at  the  beginning  of  each  time  step.  Interfaces  are  represented 
by  triangle  meshes  that  are  propagated  in  each  time  step.  The  topology  issues  of  the  interface  are  resolved  by 
the  FronTier  library  and  the  only  information  required  by  the  FronTier  library  is  the  discretized  velocity  filed  in 
the  computational  domain,  which  is  stored  in  the  center  of  each  computation  grid.  Velocity  of  each  vertex  in  the 
interface  mesh  is  the  result  of  interpolation  of  nearby  cell  center  velocities.  Then  interior  states  are  updated  by  the 
incompressible  hydro  solver. 

The  magnetic  source  term  (1(J  x  B))  is  evaluated  first.  The  discretization  of  equation  (5)  is  similar  to  that 
in  section  (3.1),  while  the  boundary  condition  is  much  simpler.  Similarly,  integrating  equation  (5)  together  with 
divergence  theorem,  we  obtain. 


which  is 


|  V-(Vip)dv=  (n  Vip-nds=  (p  -(uxB)-nJj, 

Jv  Jdv  Jdv  C 

(f)  -(u  x  B)  ■  nds. 

Jdv  C 


dip 

—  ds  = 
dv  on 


(15) 


(16) 


With  the  boundary  condition  (6),  we  can  see  that  the  integral  along  the  boundary  of  conductive  fluid  in  each  partial 
cell  is  canceled  in  both  sides  of  (16),  and  the  discretization  equation  is  greatly  simplified. 

After  solving  equation  (5),  the  gradient  of  the  electric  potential  < p  is  substituted  into  equation  (3)  and  the  current 
density  J  is  obtained.  Secondly,  we  deal  with  the  equation  (1),  without  regarding  the  divergence  constraint,  for  an 
intermediate  velocity  u*.  Employing  an  operator  splitting  technique,  we  resolve  the  advection  step 


-  =  -(u"  •  V)u".  (17) 

At 

Only  for  the  advection  step,  the  density  jump  across  the  interface  of  two  fluid  components  is  smoothed  with  a 
certain  smoothing  radius  of  computation  cells  [14].  The  advection  part,  equation  (17),  is  evaluated  explicitly, 
with  a  second  order  Godunov  type  scheme  ([10]).  For  the  diffusion  part,  we  employ  the  implicit  Crank-Nicolson 
method.  Two  fluid  components  are  solved  together,  disregarding  the  interface. 

Thirdly,  the  diffusion  step  and  the  source  term  are  resolved 


where  q  is  the  pressure  of  the  previous  time  step.  Finally,  we  perform  the  projection  step.  Applying  the  divergence 
operator  in  both  sides  of  equation 
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and  using  the  divergence  constraint  V  ■  u'!+1  =  0,  we  obtain  the  following  elliptic  equation  for  pressure 

VV+1  =  —  V  •  u*.  (20) 

At 

The  projection  step  is  an  elliptic  interface  problem  discussed  in  section  (3.1).  Two  jump  conditions  for  pressure 
are 

\P]  =  O-K 

=  0 

p  on 

where  k  is  the  curvature  and  cr  is  the  surface  tension  coefficient.  A  matrix  to  solve  such  elliptic  interface  problem 
is  set  according  to  the  algorithm  described  in  section  (3.1). 

The  pressure  is  updated  using  the  solution  of  the  projection  step: 

P  =  q  +  <f>n+1  (23) 

The  described  algorithm  achieves  the  second  order  convergence. 

The  implementation  is  carried  out  with  C++  and  MPI  for  the  communication  between  processors.  FronTier’s 
hyperbolic  solvers  demonstrate  good  scalability  on  large  machines  of  the  IBM  BlueGene  series.  The  scalability 
of  elliptic  solvers  is  determined  by  the  scalability  of  commonly  used  parallel  libraries  for  sparse  linear  system 
of  equations  (preconditioned  Krylov  subspace  iterative  solvers  of  the  PETSc  library  have  been  used  in  our  MHD 
code). 


(21) 

(22) 


4.  Verification  and  Validation 


Verification  and  validation  tests  for  the  three-dimensional  FronTier-MHD  code  have  been  performed  using 
experimental  and  theoretical  studies  of  liquid  mercury  jets  in  magnetic  fields.  Experimental  studies  of  a  mercury 
jet  entering  a  magnetic  field  with  the  magnitude  satisfying  the  hyperbolic  tangent  profile  have  been  performed 
in  [23].  An  asymptotic  theoretical  analysis  has  also  been  done  by  the  same  group.  The  experiment  setup  is  as 
follows.  A  mercury  jet  with  the  initial  diameter  of  8  mm  is  shot  horizontally  into  a  transverse  magnetic  field  with 
the  initial  velocity  of  2.1  m/s.  The  amplitude  of  the  transverse  magnetic  field  satisfies  the  following  equation 

=  -[1  -  tanh(— - — -)],  (24) 

15  max  2*  L,m 


where  zo  is  the  center  and  L,„  is  the  characteristic  length  of  the  magnetic  field.  In  our  simulations,  zo  =  1.5  cm  and 
Lm  =  0.62  cm. 

As  predicted  in  [23],  the  magnitude  of  expansion  of  the  jet  depends  on  the  z  value: 


o  ■  ,  r  cos(at) 
lisin(az  >  coshHs  t)  dt 

O  ,  '*N  r  cos(at) 

-Pco^az  )  I  coshHs  t ) dt 


(25) 


where  a  =  y/6/Wa  and /?  =  emAa/8o-.  And 

•  Na  —  (T,,B;nina/  pfUo  is  the  Stuart  number  of  the  jet,  cre  is  the  electric  conductivity  of  mercury,  a  is  the  radius 
of  the  cross-section  of  the  jet,  pf  is  the  density  of  mercury  and  cjq  is  the  main  flow  velocity  which  is  2.1  m/s. 

•  Wa  -  pfacj^/cr  is  the  Weber  number  of  the  jet,  cr  is  the  surface  tension  of  mercury. 

•  sm  =  a/Lm  and  z*  =  z/ci. 
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(a)  (b)  (c) 

Fig.  2.  FronTier-MHD  simulation  of  jet  deformation  in  magnetic  field.  Cross  sections  of  the  jet  are  shown  at  observation  points  located  at  0, 
3.5,  and  5.5  cm. 


B=1 .88T 


Fig.  3.  Mercury  jet  deformation  as  function  of  the  distance  from  the  magnetic  field  center  for  1.88  T  magnetic  field.  Results  of  simulations 
(green  dashed  line),  asymptotic  calculations  (red  dash-dotted  line),  and  experiments  (blue  dotted  line)  are  shown. 


Numerical  simulation  was  performed  using  the  magnetic  filed  strength  Bmax  of  1.41  T  and  1.88  T.  In  order 
to  save  computation  time,  simulations  were  performed  in  a  frame  moving  with  the  initial  jet  velocity.  In  the 
asymptotic  analysis  of  [23],  the  jet  was  assumed  to  extend  infinitely  and  reach  the  steady  state.  To  simulate 
similar  conditions,  initially  long  cylindrical  jet  was  moving  through  the  magnetic  field  rather  then  being  ejected 
from  the  nozzle.  Also,  the  jet  is  assumed  to  be  in  the  vacuum  while  in  the  simulation,  the  vacuum  was  substituted 
with  light  gas,  with  the  density  104  times  smaller  than  the  density  of  mercury.  With  such  a  large  density  ratio,  the 
influence  of  gas  on  the  momentum  of  the  mercury  jet  can  be  ignored.  In  order  to  obtain  accurate  profile  of  the 
electric  current  density,  the  computational  mesh  contained  approximately  20  cells  across  the  cross-section  of  the 
mercury  jet. 

Experimental  results  of  the  jet  deformation  from  [23],  results  of  asymptotic  analysis,  and  numerical  simula¬ 
tions  are  plotted  in  Figure  3  for  1.88  T  magnetic  field  and  in  Figure  4  for  1.41  T  field.  We  observe  a  very  good 
agreement  of  simulations  with  asymptotic  calculations  at  small  distances  from  the  magnetic  field  center  corre¬ 
sponding  to  smaller  jet  deformations.  The  expected  disagreement  with  experimental  results  at  small  distances  can 
be  explained  by  the  fact  that  experiments  were  carried  out  using  a  cylindrical  nozzle  located  at  z  =  0  that  reduced 
jet  deformations  compared  to  long  free  jets.  But  at  larger  distances  from  the  nozzle  corresponding  to  larger  jet 
deformations,  numerical  simulations,  theoretical  calculations,  and  experiments  are  all  in  agreement. 

We  would  like  to  comment  on  the  importance  of  maintaining  a  sharp  density  discontinuity  via  the  front  tracking 
and  embedded  boundary  methods.  Without  the  embedded  boundary  method,  the  density  ratio  that  interior  solving 
can  handle  is  limited  by  high  condition  number  of  the  corresponding  matrix  of  projection  step.  In  order  to  perform 
the  simulation  without  the  embedded  boundary  method,  we  artificially  increased  the  density  of  ambient  gas  so 
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Fig.  4.  Mercury  jet  deformation  as  function  of  the  distance  from  the  magnetic  field  center  for  1.41  T  magnetic  field.  Results  of  simulations 
(green  dashed  line),  asymptotic  calculations  (red  dash-dotted  line),  and  experiments  (blue  dotted  line)  are  shown. 


B=1 ,88T 


Fig.  5.  Degradation  of  accuracy  without  sharp  density  discontinuity.  Mercury  jet  deformation  as  function  of  the  distance  from  the  magnetic 
field  center  for  1.88  T  magnetic  field.  Results  of  untracked  simulations  (green  dashed  line),  asymptotic  calculations  (red  dash-dotted  line),  and 
experiments  (blue  dotted  line)  are  shown. 


that  the  density  ratio  dropped  to  10.  Figures  5  and  6  demonstrate  the  degradation  of  accuracy  of  simulations  if 
the  correct  density  ratio  and  sharp  density  discontinuity  are  not  resolved.  Keeping  the  discontinuity  sharp  is  even 
more  important  for  applications  involving  more  extreme  flow  regimes. 


5.  Applications 

Both  compressible  and  incompressible  fluid  FronTier-MHD  code  are  used  for  the  simulation  of  processed  rel¬ 
evant  to  energy  research  and  accelerator  applications.  Simulations  of  the  mercury  target  for  the  Muon  Accelerator 
Project  (http://map.fnal.gov)  is  among  the  most  important  applications  of  the  code.  The  target  will  contain  a  series 
of  30-cm-long  and  1-cm-diameter  mercury  jets  entering  a  strong  (~  15  Tesla)  magnetic  field  at  a  small  angle  to  the 
solenoid  axis.  When  each  jet  reaches  the  center  of  the  solenoid,  it  interacts  with  a  powerful  proton  pulse  penetrat¬ 
ing  the  jet  and  depositing  energy  of  the  order  of  100  J/g  into  mercury.  The  purpose  of  our  numerical  simulations 
is  to  evaluate  states  of  the  target  before  and  after  the  interaction  with  protons  to  optimize  the  target  design.  The 
compressible  code  deals  with  the  jet  instabilities  due  to  external  energy  deposition  and  their  partial  stabilization 
by  the  magnetic  field  [7],  The  incompressible  FronTier  MFLD  code  is  used  for  the  simulation  of  liquid  metal  jets 
in  magnetic  fields  of  different  configurations  prior  to  the  interaction  with  proton  pulses.  Other  applications  involve 
pellet  ablation  in  tokamaks  [8],  and  plasma  jet  liners  for  magneto-inertial  fusion  [9], 
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Fig.  6.  Degradation  of  accuracy  without  sharp  density  discontinuity.  Mercury  jet  deformation  as  function  of  the  distance  from  the  magnetic 
field  center  for  1.41  T  magnetic  field.  Results  of  untracked  simulations  (green  dashed  line),  asymptotic  calculations  (red  dash-dotted  line),  and 
experiments  (blue  dotted  line)  are  shown. 


6.  Conclusions 

A  numerical  algorithm  and  the  parallel  FronTier-MHD  code  for  the  study  of  three-dimensional,  incompress¬ 
ible,  multiphase  or  free  surface  MHD  flows  at  low  magnetic  Reynolds  numbers,  capable  of  handling  large  density 
ratios  across  material  interfaces,  have  been  developed.  The  algorithm  is  based  on  the  front  tracking  method  for 
material  interfaces  and  the  embedded  boundary  methods  for  elliptic  problems.  The  code  has  achieved  a  good  ac¬ 
curacy  in  verification  and  validation  tests  involving  three-dimensional  mercury  jets  in  nonuniform  magnetic  fields. 
The  FronTier-MHD  code  has  been  used  for  simulations  of  liquid  mercury  targets  for  the  proposed  muon  collider 
/  neutrino  factory,  ablation  of  pellets  in  tokamaks,  and  processes  in  hybrid  magnetoinertial  fusion. 
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